////////////////////////////////////////////////////////////////////////////////////
///////// Tempering
//////////////////////////////////////////////////////////////////////////////////


template<class Kernel, class Resample>
class Particle<Kernel, Resample, Density::Tempering>{

public:
	Particle(Kernel *K,Resample *R,Density::Tempering *D,int M,Distribution::Distribution *F){
		_p=(*D).Get_p();
		mat Xt(M,_p);//p le nombre de parametre a obtenir de Density
		double *Wt=new double[M];
		X=Xt;
		_K=K;
		_Z=0;
		_R=R;
		_F=F;
		_D=D;
		W=Wt;
		_thres=_U*M;
		_M=M;
		Normalize();
		y=X.t();
		_n=(*D).Get_n();
	//	cout << "///" << _n << "///";
	}
	void Filter()
	{
		this->Init();
		int i=1;
		double p=0; 
		//for(int i=1;i<(_n);i++)
	//	{
		double b=(*_D).Get_bn();
	//	b+=(double)(1/(double)(_P-1));
		double temp=StepLength(b);
		b+=temp;
		cout << "b: " << b << " p";
		_D->Set_Phi(b);
		y=X.t();
		_D->Weight(y,W,i);//le p designera comment on avance 
		cout << " ESS: " << Ess_W(W,1);
		X=y.t();
		while(p<1){
		//	cout << i;
			cout << "step:" << i << "\n";
			this->Step(i);
			p=_D->Get_bn();
		}
		
	}
	void Init()
	{
		//pior
		for(int i=0;i<_M;i++)
		{
			W[i]=1;
			//mat A=p.r();
			X(i,span::all)=(_F->r(i)).t();
		}
	//	cout << "Hello";	
	}
	void Step(int i)
	{
		/*double Ess=Ess_W(W);
		//cout << "Ess:" << Ess << "\n";
		double b=_D->Get_bn();
	//	b+=(double)(1/(double)(_P-1));
		double temp=StepLength(b);
		b+=temp;
		cout << "b: " << b << "p";
		_D->Set_Phi(b);
		//cout << "\\";
		*/
		mat v=growingvect(_M);
		double *w=new double[_M];
		for(int j=0;j<_M;j++){
			 w[j]=W[j];
		//	cout << w[j] << "\n";
		}
		(*_R)(&_M,w,v);
		//cout << v << "\n";
		delete[] w;
		Arangemat(X,v);
		Normalize();//poids fixé a 1
	//	mat C=cov(X);
		//cout << C;
		cout << "sum: "  << sum(X,0) << "\n";
		_K->Set_s(X);
		//cout << "before:" << X.row(1);
//		for(int k=0;k<3;k++)
//		{	
			_K->Move(&(X),i);
//		}
		//cout << "after:" << X.row(1);
		double b=_D->Get_bn();
	//	b+=(double)(1/(double)(_P-1));
		double temp=StepLength(b);
		b+=temp;
		cout << "b: " << b << "p";
		_D->Set_Phi(b);

		y=X.t();
		//cout<< "before:" << W[0] ;
		//cout<< "after:" << W[0] ;
		//dans le cas de tempering peut etre determiner
		_D->Weight(y,W,i);//le p designera comment on avance 
	//	cout << "after:" << y.col(1);
		X=y.t();
		cout << " ESS: " << Ess_W(W,1);	
	}

	double Ess_W(double *w,int tf)//renormalisation???
	{
		double *foo= new double[_M];
		double sum=0;

		for(int i=0;i<_M;i++)
		{
			sum+=w[i];
		}
		if(tf){
			cout << "SUM EV : " << sum;	
			_Z+=log(sum/_M);
		}

		for(int j=0;j<_M;j++)
		{
			foo[j]=(w[j]/sum);
		}

		sum=0;
		for(int i=0;i<_M;i++)
		{
			sum+=foo[i]*foo[i];
		}
		delete[] foo;

		return 1/sum;		
	}
	double StepLength(double phi)
	{
		double u=1.05-phi;
		double l=0;
		double thres=0.001;
		double e=2;
		double alpha=0.05;
		double eta=0;	
		while((e>thres)&(l<1-phi))
		{
			mat y=X.t();
			eta=Ess_W(_D->Eta_star(alpha,W,y),0)/_M;
			if(eta>_U)
			{
				l=alpha;
				alpha=(double)(alpha+u)/2;
			}else{
				u=alpha;
				alpha=(double)(alpha+l)/2;

			}	
			e=abs<double>(u-l);
			cout << "\\\\\\ "<< eta <<" \\\\\\\\\n"; 
		}
		return min2<double>(alpha,1-phi);
	}

	
	mat Get_theta(void){return X;}
	void Normalize(void){for(int i=0;i<_M;i++){W[i]=1;}}
	mat   Get_W(void){
		double sum=0;
		mat C(_M,1);
		for(int i=0;i<_M;i++)
		{
			sum+=W[i];
		}	
		for(int i=0;i<_M;i++)
		{
			C(i,0)=(W[i]/sum)*_M;
		}
//		cout << C;
		return C;
	}
	double Get_Z(void){
		return _Z;}

private:
	mat X;//value particle
	double *W;//weigths
	double Ess;
	int _M;
	double _thres;
	double _Z;
	Distribution::Distribution *_F;
	Kernel *_K;
	Resample *_R;
	Density::Tempering *_D;	
	mat y;
	int _n;
	int _p;

};

